
clear all;

A = { "gsl_ext_conjgrad_fr_mt", "gsl_ext_conjgrad_pr_mt", ...
      "gsl_steepest_descent", "gsl_ext_bfgs_mt", ...
      "gsl_ext_lrwwsimplex", "gsl_ext_mnewton", ...
      "gsl_conjugate_pr", "gsl_vector_bfgs" };

NUM_ALGORITHMS = length(A);

MIN_ALPHA_EXPONENT = -6;
MAX_ALPHA_EXPONENT = 0;
NUM_ALPHA_VALUES = 25;
ALPHA_VALUES = logspace(MIN_ALPHA_EXPONENT, MAX_ALPHA_EXPONENT, NUM_ALPHA_VALUES);

MAX_NUM_ITER = 1000;

stopcrit_params = struct;
stopcrit_params.xmin = [ 0, 0 ]';

for ai = 1:NUM_ALGORITHMS
  R = NaN(NUM_ALPHA_VALUES, 1);
  for alphai = 1:NUM_ALPHA_VALUES
    alpha = ALPHA_VALUES(alphai);

    f = [ '(', num2str(alpha, '%.20f'), '*x)^2+y^2' ];
    x0 = [ 3 / alpha, 2.1 ];
    stopcrit_params.eps = 1 / alpha * 1e-8;

    algoparams = struct;
    algoparams.initsimplexsize = [ 1 / alpha, 1 ]';

    results = GSLpp_minimize(A{ai}, ...
              algoparams, struct, f, struct, ...
              struct('f', 'sym', 'g', 'sym', 'H', 'sym'), ...
              "xdisttomin", stopcrit_params, ...
              x0, 0, false);

    if results.converged == true && results.numiter < MAX_NUM_ITER
      R(alphai) = results.numiter;
    endif
  endfor

  minalphai = NUM_ALPHA_VALUES;
  maxnumiter = 0;
  for alphai = 1:NUM_ALPHA_VALUES
    if ~isnan(R(alphai))
      minalphai = min(alphai, minalphai);
      maxnumiter = max(maxnumiter, R(alphai));
    endif
  endfor

  GSLpp_newplot;
  semilogx(ALPHA_VALUES, R, 'b', 'linewidth', 6);
  axis([ ALPHA_VALUES(minalphai), ALPHA_VALUES(NUM_ALPHA_VALUES), 1, maxnumiter ]);

  title(GSLpp_title(A{ai}));
  xlabel('\alpha');
  ylabel('Number of iterations');

  print([ "param_scaling_test_", A{ai}, ".ps" ], "-landscape", "-dashed", "-FArial:23");
  close;
endfor
